home *** CD-ROM | disk | FTP | other *** search
/ FM Towns: Free Software Collection 6 / FM Towns Free Software Collection 6.iso / ms_dos / happypas / pai.pas < prev   
Pascal/Delphi Source File  |  1993-07-08  |  5KB  |  160 lines
  1. {*********************************************************************
  2.  *        円周率の計算                                               *
  3.  *             HAPPyのサンプルプログラム (TOWNSフリコレ6版)          *
  4.  *             Auther  H.Asano                                       *
  5.  *********************************************************************}
  6.  
  7. { 定数AREA-7桁分の円周率を求めます。
  8.   しかし、不思議なことに、小数点以下11桁目の答えが違うのです。
  9.   正解は8なのに、このプログラムでは9となります。
  10.   1000桁まで計算すると、761,762桁目も違うことが判明しました。
  11.   どなたか、このプログラムの誤りを指摘して下さい  }
  12.  
  13. { マーチンの公式 π/4 := 4*arctan(1/5) - arctan(1/239) を 使う。
  14.     π   := 16*arctan(1/5) - 4*arctan(1/239)
  15.     arctan(x) := x/1 - x*x*x/3 + x*x*x*x*x/5 - ･･･
  16.     x/1の分子分母にxをかけて整理すると
  17.     16*arctan(1/5)   := 16*5/25/1 - (term)/25/3 + (term)/25/5 - ････
  18.                       = 80  /25/1 - ･･･
  19.     -4*arctan(1/239) := -4*239/57121/1 + (term)/57121/3 - (term)/57121/5 + ･･･
  20.                       = -956  /57121/1 + ･･･
  21.     ここで、termは、前項までのx指数計算の結果を格納している。
  22.     たとえば第1項の80/25が次の項のtermである。第2項は 80/25/25 / 3 となる。
  23. }
  24.  
  25. program Pai(output) ;
  26.  
  27. const AREA = 187 ;                         { 180桁まで求める }
  28.      {***** このプログラムをHAPPyのスピードアップ測定に使ってきました。
  29.             使用マシンは、FM TOWNS 40Hです。
  30.              以下は､スピードアップの履歴です                    *****}
  31.                   { 1992.4.13                      　1分55秒 }
  32.                   { 1992.4.18                      　1分49秒 }
  33.                   { 1992.4.18 通訳系の最適化により   1分44秒 }
  34.                   { 1992.4.22 通訳系の書き直しにより 1分14秒 }
  35.                   { 1992.4.30 直訳系と通訳系の分離   1分10秒 }
  36.                   { 1992.5.01 通訳系の記憶装置型変更 1分00秒 }
  37.                   { 1992.5.02 通訳系の改良           0分58秒 }
  38.                   { 1992.5.02 program の改良         0分56秒 }
  39.                   { 1992.11.03 通訳系の改良          0分54秒 }
  40.                   { 1992.11.26 通訳系の改良          0分48秒 }
  41.                   { 1993.02.17 programの改良         0分38秒 }
  42.  
  43. type range   = 0..AREA ;
  44.      wkplace = array[range] of integer ; { 添字2と3の間が小数点 }
  45.                                 
  46. var  pai          : wkplace ; { 円周率の格納エリア }
  47.      work,term    : wkplace ; { 計算上のワーク     }
  48.      i            : range   ;
  49.      p            : range   ; { 収束位置 }
  50.  
  51. {********** 割り算処理 **********}
  52. {  wk := term div jyosu を 計算する。
  53.    pは、termi配列の0でない最初の位置を示す。これは収束判定に使う }
  54. procedure Divide(var wk : wkplace ;  jyosu  : integer) ;
  55.   var i    : range   ;
  56.       hi,w : integer ;
  57. begin
  58.                   { 有効数字が出てくる位置を求める  }
  59.   while (term[p] = 0) and (p < AREA) do p := succ(p) ;
  60.  
  61.   hi := 0 ;
  62.   for i := p  to  AREA do 
  63.   begin
  64.     hi := hi*10 + term[i] ;
  65.     w  := hi div jyosu    ;
  66.     if w >= 1 then
  67.     begin wk[i] := w ;
  68.           hi    := hi mod jyosu
  69.     end 
  70.     else  wk[i] := 0
  71.   end
  72. end ;
  73.  
  74. {********** 加減算 **********}
  75. {  pai := pai + work  ････ AddSubFlag = true  の時 
  76.    pai := pai - work  ････ AddSubFlag = false の時 }
  77. procedure  AddSub(AddSubFlag  : Boolean) ;
  78.   var carry : Boolean ;       { 次の桁に影響する時 true }
  79.       wk    : integer ;
  80.       i     : range   ;
  81. begin
  82.   carry := false ;
  83.  
  84.   if AddSubFlag then
  85.   begin                             { 加算 }
  86.     for i:=AREA downto p do
  87.     begin
  88.       wk := pai[i] + work[i] + ord(carry) ;
  89.       carry := wk >= 10 ;
  90.       if carry then pai[i] := wk - 10
  91.                else pai[i] := wk 
  92.     end
  93.   end
  94.  
  95.   else
  96.   begin                             { 減算 }
  97.     for i:=AREA downto p do
  98.     begin 
  99.       wk := pai[i] - (work[i] + ord(carry)) ;
  100.       carry := wk < 0 ;
  101.       if carry then pai[i] := wk + 10
  102.                else pai[i] := wk
  103.     end    
  104.   end
  105. end ;
  106.  
  107. {********** 結果の印字処理 **********}
  108. procedure Print  ;
  109.   var  i : range ;
  110.        j : integer ;
  111. begin
  112.   j := 0;
  113.   write('円周率= ',pai[2]:1,'.');
  114.   for i:=3 to AREA-5 do 
  115.   begin
  116.     j := j + 1 ;
  117.     write(pai[i]:1) ;
  118.     if j mod 5 = 0 then  write(' ') ; {  5桁ごとに空白 }
  119.     if j=50 then                      { 50桁ごとに改行 }
  120.     begin
  121.       writeln ;
  122.       write(' ':10) ;
  123.       j:=0 
  124.     end 
  125.   end ;
  126.   writeln('･･･')
  127. end ;
  128.  
  129. {********** 計算処理 **********}
  130. procedure calc(j1 : integer ; Plus : Boolean) ;
  131.   var j2    : integer ;
  132.       add   : Boolean ;
  133. begin
  134.   p   := 0 ;
  135.   j2  := 1 ;
  136.   add := plus ;
  137.   while p < AREA do
  138.   begin
  139.     Divide(term,j1) ;        { term := term / j1   }
  140.     Divide(work,j2) ;        { work := term / j2   }
  141.     AddSub(add) ;            { pai  := pai +- work }
  142.     j2  := j2 + 2 ;
  143.     add := not add
  144.   end
  145. end ;
  146.  
  147. {********** メイン処理 **********}
  148. begin {main}
  149.   for i:=0 to AREA do pai[i] := 0 ;
  150.  
  151.   term := pai ; term[1] := 8 ;                      { 初期値  80.0 }
  152.   calc(25{5*5},true{+-+-････}) ;        { マーチンの公式 第1項計算 }
  153.  
  154.   term[0] := 9 ; term[1] := 5 ; term[2] := 6 ;
  155.   for i:=3 to AREA do term[i] := 0 ;                { 初期値 956.0 }
  156.   calc(57121{239*239},false{-+-+･･･}) ; { マーチンの公式 第2項計算 }
  157.  
  158.   Print
  159. end.
  160.